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Abstract 

Background: Most of the current biophysical models designed to address the large-scale distribution of malaria 
assume that transmission of the disease is independent of the vector involved. Another common assumption in these 
type of model is that the mortality rate of mosquitoes is constant over their life span and that their dispersion is 
negligible. Mosquito models are important in the prediction of malaria and hence there is a need for a realistic 
representation of the vectors involved. 

Results: We construct a biophysical model including two competing species, Anopheles gambiae s.s. and Anopheles 
arabiensis. Sensitivity analysis highlight the importance of relative humidity and mosquito size, the initial conditions 
and dispersion, and a rarely used parameter, the probability of finding blood. We also show that the assumption of 
exponential mortality of adult mosquitoes does not match the observed data, and suggest that an age dimension can 
overcome this problem. 

Conclusions: This study highlights some of the assumptions commonly used when constructing mosquito-malaria 
models and presents a realistic model of An. gambiae s.s. and An. arabiensis and their interaction. This new mosquito 
model, OMaWa, can improve our understanding of the dynamics of these vectors, which in turn can be used to 
understand the dynamics of malaria. 
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Background 

This is the first of two papers describing a dynamic 
model (Open Malaria Warning; OMaWa) of Anopheles 
arabiensis and Anopheles gambiae s.s. Our aims in 
this article are 1) to formulate recent research on the 
Anopheles gambiae complex in a mathematical frame- 
work, and 2) to show how the new formulations influence 
the dynamics of malaria and mosquito populations. 

In this paper, we describe a model of the dynamics 
of the two species and then show how parameters can 
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influence the success of the two species, and how temper- 
ature, humidity and mosquito size can influence malaria 
transmission. 

Climate and malaria 

Most of the 149-274 million cases and 537,000-907,000 
deaths from malaria occur in sub-Saharan Africa [1,2]. 
Climate has been one of the main drivers of this dis- 
ease [3], governing the spatial extent and year-to-year 
variations. The pathway from climate to malaria goes 
through the parasite and the mosquito. Although it is well 
established [4] how parasite development is influenced by 
temperature [5], the vectors response to weather and cli- 
mate is more complex. Mosquito density depends not only 
on temperature but also on the abundance of breeding 
sites (rainfall and evaporation) [6], desiccation (humidity) 
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[7], and competition between mosquitoes [8]. In the past 
20 years, a shift in the distribution of An, arabiensis and 
An. gambiae s.s, has been observed in Kenya [9], show- 
ing that the species composition is not static over time. 
In the context of climate change [10], variability in vec- 
tor populations is a factor that has not been considered 
so far. 

Malaria and mosquito models 

At the turn of the 20th century the work of several 
researchers, including Battista Grassi and Ronald Ross, 
resulted in the discovery that mosquitoes of the Anopheles 
genus transmit malaria [11,12]. Over the next 20 years, 
Ross, and later Lotka and Waite, developed mathemati- 
cal models that became central in malaria control [13-19]. 
In the 1950s, George MacDonald refined these models 
and showed that DDT could be used to interrupt malaria 
transmission [20]. Since then, several modelers have fol- 
lowed in the footprints of Ross, Lotka, and MacDonald 
[21-30]. Some have designed models to show how tem- 
perature alone influences malaria transmission [31], while 
others have focused on the theoretical eff^ect of bed 
nets [32], multiple interventions [33] or climate change 
[34-36]. There is also a growing number of models that 
address the dynamics of immunity within individuals [37] 
and in communities [21,38]. 

In 2011, The malERA Consultative Group on Modeling 
[39] provided a review of the current state of mathematical 
models and pointed to the importance of good mosquito 
models for assessing the impact of climate change on 
malaria. 

Many traditional models rely on a threshold principle. 
The idea has been to find thresholds for longevity, num- 
ber of bites or days to recovery that must be reduced 
to interrupt the transmission. With increased computa- 
tional power it is now possible to make more complex 
models and hence explore a wider range for the dynam- 
ics of malaria and mosquito survival. By integrating the 
knowledge from simpler models into a complex system, 
it is possible to test if the assumptions are true over 
a wider geographical range. In addition, these complex 
models can make quantitative predictions about strategies 
for control [40]. 

Model summary and motivation 

A model is mental copy that describes one possible rep- 
resentation of a system. We present an alternative for- 
mulation of the dynamics of An, gambiae s.s. and An. 
arabiensis. The model is a system of ordinary differen- 
tial equations (ODEs) with three compartments: eggs, 
first to fourth instar larvae, and pupae; an age-structured 
formulation of adult mosquitoes; and size prediction for 
adult mosquitoes (measured as wing length in mm). 
This can be considered the skeleton of the model. As 



demonstrated later, the model structure can be simplified 
when mosquito size can be neglected or when we assume 
no births. The model can be run with a spatial structure in 
which we include or exclude mosquito dispersion, or as an 
idealized model in which the model is evaluated at a single 
point. 

The ODEs parametrize daily mortality rates, which are 
size-dependent for adult mosquitoes; development rates 
in the aquatic stages; biting rates; fecundity; the proba- 
bility of finding a blood meal; and mortality related to 
flushing of eggs, larva and pupa out of oviposition sites. 
These parametrization schemes are driven by air tem- 
perature, relative humidity, relative soil moisture, water 
temperature, and runoff. As already mentioned, the model 
can be applied in a spatial domain. In this case, tem- 
perature and other environmental data are taken from a 
regional climate model, the Weather Research and Fore- 
casting Model (WRF) [41]. In the examples shown later, 
we run the model at a resolution of approximately 50 
km and a temporal resolution of 5-20 years in steps of 
3 h. In addition to weather data, human [42] and cattle 
[43,44] densities are introduced to estimate the probability 
of feeding. 

At this spatial resolution, the model should potentially 
be able to define larger foci of mosquito productivity, 
while the ability to identify hotspots will be limited [45] . 
However, 50 km is the standard for regional climate mod- 
els addressing long-term changes in climate [46]. In addi- 
tion, the true accuracy of historical cattle and human 
population density estimates for Africa in general is not 
likely to be greater than 50 km. 

The mosquito model described here is designed to cap- 
ture the spatial distribution and the time-dependent den- 
sity of An. gambiae s.s. and An. arabiensis. If the model 
can capture the current distribution and density of the 
two species and how they are related to malaria, a future 
version of this model, including infections, could be used 
to explore the long-term impact of current interventions 
under a changing climate. To have confidence that the 
model has these abilities, several aspects not considered 
here should be evaluated (papers under preparation). In 
addition, if malaria modelers move towards the ensemble 
thinking widely adopted in the climate community, this 
model could be one representation of historical and future 
changes for malaria. The aim of such an ensemble would 
be to deal with uncertainties in the system. Ultimately, 
the goal would be to produce policy-relevant information 
including uncertainty. 

We have chosen to represent the non-exponential mor- 
tality of An. gambiae s.s. and An. arabiensis as observed 
in laboratory settings [47], semi-field conditions [48], and 
in the field [49]. A common assumption is that in the 
field, mortality rates are constant with age because of 
predation [31]. To date, few studies have confirmed this, 
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while there is field-based evidence of age-dependent Ae. 
aegypti mortality [49], which has implications for malaria 
transmission [50]. In the model, we also describe how 
mosquito size changes over the season. This might seem 
to be an overcomplication of the model. The motivation, 
however, is that we have observed substantial improve- 
ments for arid regions such the Sahel when we included 
mosquito size prediction. Fouet et al. reported that 
mosquito size is an important adaptation strategy in arid 
environments [51]. 

We do not claim that the additional complexity adds 
any value. Stating this before the model has been fully 
evaluated and compared to simpler models would be 
dangerous. The model is thus one possible way of describ- 
ing the dynamics of An. gambiae s,s, and An, ara- 
biensis. It is under continuous development, and we 
expect to add and alter components as new data become 
available. 

To highlight some of the components that contribute 
to the dynamics of An. gambiae s.s. and An. arabiensis in 
the model, five sensitivity experiments focus on the effect 
of temperature, relative humidity and mosquito size on 
malaria transmission. We also show how An. gambiae s.s. 
and An. arabiensis respond to changes in the probability 
of finding blood, carrying capacity, initial conditions, and 
dispersion. 

Material and methods: model description 

Summary of the model 

Figure 1 provides an overview of the model. In the follow- 
ing sections we present the ideas behind the model and 
its general structure, how a climate model is used to drive 
the mosquito model, and the parametrization schemes 
used in the model. It should be possible to read each part 
independently; for example, data from a climate model 
can be used to drive any malaria model; the parametriza- 
tion scheme can be used in any malaria model; and the 
malaria model described here can be used with differ- 
ent parametrization schemes, with or without data from a 
climate model 

As mentioned above, the model comprises a system of 
ODEs for eggs, first to fourth instar larvae, and pupae; 
an age-structured formulation for adult mosquitoes; and 
size prediction for adult mosquitoes (measured as wing 
length in mm). The first limitation in the aquatic stage is 
the availability of ovipositing sites, which is parametrized 
in terms of relative soil moisture and the potential for 
puddle formation in a specific location. Once oviposit- 
ing sites have been formed, adult female mosquitoes are 
allowed to deposit eggs until the site is full, defined 
as the biomass relative to the carrying capacity for the 
location. To account for density-dependent mortality, 
first instar larvae can be preyed on by fourth instar 
larvae [52], and an extra density-dependent mortality 



term is added to account for prey-independent mortal- 
ity [53]. The numbers of eggs, larvae and pupae are 
reduced when the precipitation rate exceeds the infil- 
tration rate. The larval density in the aquatic habi- 
tat influences the size of adult mosquitoes [53]. We 
account for this by predicting mosquito size at emer- 
gence as a function of larval density. In addition to 
temperature and relative humidity [47], mosquito size 
influences the daily adult survival probability [7,51,54,55] 
([56], Aedes aegypti). We therefore describe an adult 
survival model that takes temperature, relative humid- 
ity and mosquito size into consideration. In addition, 
adult mortality and fecundity can increase if there are 
no or few sources of blood. This follows the idea that a 
mosquito living in an environment where much energy 
has to be used to find blood will do this at the cost of 
survival. 

We adopt these general ideas for two species, An. 
gambiae s.s. and An. arabiensis. It should be noted that 
we have less confidence in the model for the An. gambiae 
s.s. M form, since aestivation (as documented by Lehmann 
et al. [57] and Adamou et al. [58]) is not included. 
In addition, there are some indications that the M 
form breeds in larger pools [59] and hence the pud- 
dle parametrization might have limited validity for this 
form. 

In addition to time, the model can include two (three, 
since space is two-dimensional) additional dimensions, 
namely age and space. The space dimension allows dis- 
persion of mosquitoes, meaning that (re)establishment 
through migration to areas that were previously free of 
An. gambiae s.l. is possible. The gradual invasion of Brazil 
by An. arabiensis in the 1930s [60] is one example of 
dispersion. 

The ODEs were solved using the ODE solver Isoda 
[61-63]. The relative and absolute error tolerances were 
not modified from the original Isoda implementation 
(le~^). The model can be run either as a spatial model 
(with or without mosquito dispersion) or evaluated at a 
single point at which movement is neglected. A detailed 
overview of the possible model parameters can be found 
in Table 1. 

Differential equations for the aquatic compartment 

The aquatic compartment consists of six stages: eggs 
(£), four larval stages (Li,L2>^3>^4)> and pupae {P). 
Transitions between the different compartments can be 
expressed in terms of delayed equations. To simplify the 
solution and avoid numerical instabilities, we approximate 
the model as ODEs [21]. Lunde et al. reported on the 
errors introduced by this approximation [64]. 

New eggs added to the population depend on the num- 
ber of adult mosquitoes (m), the size of adult mosquitoes 
{i^size)y the inverse length of the gonotrophic cycle (G(r)), 
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Figure 1 Overview of tiie mosquito model. A (regional) climate model is used to force the mosquito model. In addition, static and semi-static 
fields are used as part of the parametrization schemes. Human and bovine densities limit the availability of blood meals. 



how much water is available {SMr, dimensionless) and the 
larval biomass already present in puddles (Bi): 

^ = 6 (m, msi,e.) • G{T) • SMr • (^1 " ^ j (i) 

-{PN,E(T)+Pi,E + rE(T))'E, 

where ^(m.msizen) represents potential new eggs from 
each age group, G(T) is either constant or dependent on 
temperature T, SMr is a function of the relative soil mois- 
ture and the potential puddle formation area, K is the 
maximum larval biomass a grid cell can hold, Pn,e(T) is 
natural mortality rate for eggs [Eqs. (16) and (18)], Pj^e is 
the induced mortality rate for eggs (not specified) and te 
is the inverse of development time from eggs to first instar 
larvae. 

The term 1 — Bl/K is used as a scaling factor to mod- 
ify the growth rate. When the population is low compared 
to the breeding sites available, its growth is high. As the 
population grows, there is more competition for food, 
predators become more abundant, and the growth slows. 
In the egg compartment this represents the idea that the 
mosquitoes will lay fewer eggs when breeding sites are 
already occupied [65]. 

First instar larvae (Li) are added as eggs develop into lar- 
vae. Additional mortality is added in the transition stage 
in relation to how much biomass there already is in a given 
location [53]. This approximation of increased (density- 
dependent) mortality arises because of competition and 
predators; if a puddle already is full, the number of eggs 



developing to first instar larvae is reduced, whereas if a 
puddle is empty (1 —Bl/K = 1), no extra mortality occurs. 
Similar terms could have been added to the second, third 
and fourth instar larvae, but we assume that earlier life 
stages will be affected more by density-dependent compe- 
tition and predation. 

Shoukry looked at how fourth instar larvae of An. 
pharoensis prey on first instar larvae during a 24-h exper- 
iment [52]. Using these data, we add additional mortality 
for first instar larvae according to the density of fourth 
over first instar larvae. The constant Cpred is tunable to 
both limit the predation on L\ and make it more spe- 
cific to species in the future. At most temperatures, this 
constant does not influence the density of mosquitoes 
(Additional file 1). 

The number of first instar larva is given by: 



^ = r.(r).£.fl-^ 



0.4465 



^-{liN.L{T) + lii,L + TL,{T)) 
2.989T ' ^Pred- 



(2) 



Second (L2), third (L3) and fourth instar larvae (L4) and 
pupae [P) are controlled by the development rate r and 
mortality fi: 

^ = XL,{T)-Lr - {Pn,l(T) + Pi,L + ri,(r)) .I2 (3) 
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Table 1 Model parameters 

Variable Description 



Equation(s)/reference 



Tjndoor 



rrin 
PiB) 

water 
Tsoil 

^N,L (Twater) 
T^gamb 

TE 

rp 

forob 
fgamb 

Ln 

F arab 
F gamb 
Sf 

Frm 
E 

G(T) 
t 

Bl 

SMr 



K 

/-I 

Ls 
La 
P 



Indoor temperature 

Near surface 
temperature (2 m) 

Potential number of 
new eggs 

Number of mosquitoes in 
each age group 

Daily probability of getting 
a blood meal 

Water temperature 

0-10 cm soil temperature 

Natural mortiality rate, 
eggs, larva, and pupa 

An. gambiaes.s. develop- 
ment rate, aquatic stages 

An. orobiensis develop- 
ment rate, aquatic stages 

An. gambiae s.l. develop- 
ment rate, eggs 

An. gambiae s.l. develop- 
ment rate, instar 1-4 

An. gambiae s.l. develop- 
ment rate, pupa 

Aquatic development rate 
modification An. arabiensis 

Aquatic development rate 
modification An. gambiae 

5.5. 

Number of larvae 

Mortality rate modification 

Mortality rate modification 

scaling factor for wind 
dispersion 

Flight range 

Number of eggs 

Biting rate/gonotrophic 
cycle 

time 

Larva biomass 

Induced mortality 

in aquatic and adult stages 

Dimensionless time vary- 
ing water constant, or rate 
at which ovipositing sites 
are found 

Carrying capacity 

Number of 1^^ instar larva 

Number of 2"^ instar larva 

Number of 3^^ instar larva 

Number of 4^^^ instar larva 

Number of pupa 



36 

25, 26, 30, 36 
13 



41 

14,16,18 
[91-94] 
14,1,2, 3,4, 5,6 

20 

22 

[97] 1 
[97] 2, 3,4,5 

[97] 6 
[8] 



21, 19 
[72] 17 
[72] 15 
39 

41 
1 

26 



1,2, 3,4, 5,6,7, 



24 



24 
2 
3 
4 
5 



Table 1 Model parameters (Continued) 

Cpred Predation constant. 

Currently set to 0 



' gonot 

Dd 

Tc 

h 
n 

msize 



ER 

D 
LT 

K 

HBI 



' mod 

Pbovine/cattle 
Phuman 



part of gonotrophic cycle 
formulation 

Degree days 

Critical temperature 

Adult mortality related to 
feeding 

Number of humans 
flux 

Dimension in age grid 

Size of newly emerged 
mosquitoes 

Size of mosquitoes in age 
group n 

Prediction of larva size 
Size constant 
Size constant 

Potential river length in km 

Equally spaced river 
dataset resolution in 
degrees 

Earth radius in 
km (6371.22) 

latitude in radians 

Diffusion coefficient 

Local time 

Diurnal modification for 
transport of mosquitoes 

Human blood index 

Size dependent mortality 

Natural mortality of adult 
mosquitoes 

Survival curve for adult 
mosquitoes 

Shape parameter for adult 
survival 

Sub-function for 
equation 33 

Probability of finding cattle 

Probability of finding 
humans 



Description of model components. 



2 
26 

[108], 26 
26 
42 

[42] 
39 

9 

12 

10 
[22] 
[22] 
23 
23 

23 

23 
39 
37 
37 

41,42 
28 
32, 7, 8 

35,31 

3330 

34 

41 
41 



8L3 

8t 
8P 



= XL,{T) .12 - {fiN,L{T) + fii,L + TL,{T)) • I3 (4) 

= XL,{T) .I3 - {Pn,l{T) + Pi^L + rL,{T)) -U (5) 
= XL,{T) • U - {PnAT) + Pi,p + rp{T)) . (6) 
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where P is the daily mortality rate, with the first subscript 
denoting natural (N) or induced (I) mortality and the sec- 
ond subscript denoting the aquatic stage. The subscript 
for the development rate, r, corresponds to the aquatic 
stage. The parametrization schemes and data sources used 
to estimate the rate at which eggs are laid {G{T) and 6), 
mortality (^) and the development rate (r) are discussed 
later. 

Differential equations for adult mosquitoes 

The life history and mortality rate vary over the lifespan of 
a mosquito population. We formulated a model to account 
for this variation. Adult mosquitoes are denoted by m^, 
where n indicates the age group; fz = 1 is the youngest 
group and n = 9 refers to the oldest mosquitoes. The 
age groups in the model are mi =[0,1], m2 = (2,4], 
ms = (5,8], m4 = (9,13], ms = (14,19], me = (20,26], 
my = (27,34], ms = (35,43] and m9 = (44, oo] days, 
with ageing coefficients an of 1.000, 0.500, 0.333, 0.250, 
0.200, 0.167, 0.143, 0.125 and 0.067 for n = 1,2,..., 9, 
respectively. Mosquito ageing is represented by 
where n denotes the age group. Ageing is time-invariant 
and is thus not related to the number of gonotrophic 
cycles. 

Although there is no ageing from age group 9, the 
term ^9 is included to limit the concentration of old 
mosquitoes. This is a user-specified variable and in the 
model results shown here we set this to Ys^^y~^ 
arabiensis and An, gambiae s.s.; this value should be set 
to ensure that mosquito populations can survive during 
dry periods [66,67], but still hinder accumulation of old 
mosquitoes. This can be particularly useful if the mortal- 
ity model described later is replaced with a model in which 
mortality is independent of age. 

When m is written with subscripts i and j in addition 
to Hy this denotes inclusion of mosquitoes from neigh- 
boring areas. For example, subscript i — 1 indicates that 
mosquitoes to the west of the point of interest are interact- 
ing with the point of interest. The formulation presented 
here includes movement of mosquitoes, and where appro- 
priate we denote mosquitoes by mn,i,j^ 

Again, p denotes mortality, with the first subscript 
denoting natural (N) or induced (I) mortality and the 
second subscript denoting the age group (m«) of the 
mosquitoes. O represents the mosquito flux (transport) 
and subscripts i and j define which boundaries are eval- 
uated. This is discussed in the section "Movement of 
mosquitoes". 

The number of adult mosquitoes of a specific age 
in a grid point is controlled by new mosquitoes from 
myi-i, as well as the flux to and from the point of inter- 
est {^]=-i Y}j=-i ^ij^n^i^y natural mortality PN.mn^ 
induced mortality ageing to m^+i, and mortality 



due to lack of food {P(B)), Parametrization schemes 
related to mortality are discussed later. 

This results in the following equation for the first age 
group: 

i=-lj=-l V ) 

- {PN,mi + Pl,mi + ^1) • ^1. 

The equations for age groups n =[2,9] are 
8mn ^ ^ 

i=-lj=-l 

Differential equations predicting mosquito size 

Mosquito size (msize) is important for the efficiency of 
mosquito multiplication. There are also some indications 
that increased body size is a strategy for survival in arid 
environments [7]. In general, high larval density leads to 
a smaller body size as adults, and vice verse [68]. Where 
only one species is competing for a resource, such as in 
a small puddle, mosquito size, and hence the number of 
eggs laid by each mosquito, will be of less importance. 
If two species are competing for the same resource (e.g. 
An. arabiensis and gambiae s.s,), the trade off between 
development time and size can be important in competi- 
tion for breeding sites. An, gambiae s,s, generally develop 
faster than An, arabiensis, but end up with a smaller 
body size. An, arabiensis spends more time in the aquatic 
stages and develops larger bodies, and can thus produce 
more eggs. Since our model includes competition between 
those species, we describe mosquito size as a function 
of competition for breeding sites. In theory this should 
improve our ability to separate geographical and seasonal 
distributions of An, arabiensis and An, gambiae s,s. 

Since the size of An, arabiensis and An, gambiae s,s, 
stabilizes after approximately 4 days [7] and ovoposition 
does not start before this, it is not necessary to differen- 
tiate the maximum and minimum size depending on age 
to mimic changes in the number of eggs per mosquito 
with age. However, this may be required if mortality based 
on desiccation [7,69] is used. Although mosquito size at 
a given time can be approximated using finite differences, 
we develop a different approach that is more efficient in 
terms of computational time in our model framework. 
Mosquito size for the first age group depends on larval 
size. Since the pupation time is short, this assumption 
is justified, although it might introduce minor errors. In 
a future version of the model, we plan to predict lar- 
val size dynamically. The limitations set on mosquito size 
(described in "Parametrization schemes in the aquatic 
stages") in this model might lead to An, arabiensis that 
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are slightly too small compared size in the field study of 
Ye-Ebiyo et al. [70], but the size is in line with studies 
by Huestis et al. [71] and Kirby et al. [72]. Kirby et al. 
also noted that mixed populations of An. arabiensis and 
An, gambiae s,s, had a negative effect on mosquito size 
at some temperatures. This mechanism is not included 
in the current work. However, the most important aspect 
of modelling of mosquito size is to capture seasonal and 
spatial variations. 

For size prediction we use the symbol nisizeny where n is 
the age group as described above. 

The size (wing length in mm) of newly emerged 
mosquitoes is approximated according to the linear 
relationship 

^sizce ~ 1«25 + 5 • Lsizet (9) 

where larva size Lsize (in mg) is approximated as: 

(10) 



The constants aspp and bspp are 0.45 and 0.12 for 
An, arabiensis and 0.383 and 0.147 for An, gambiae s,s,, 
respectively [22]. 

The size of mosquitoes in the first age group at any time 
is given by 



8m 



sizei 



8t 



= mm max 



J I ■■'■sizee . 

'log \ ) -msizei- 

^sizei 



(11) 



Therefore, the size of newly emerged mosquitoes 
{y^sizei ) depends on the number of newly emerged pupae 
and the relative density of larva at the breeding site. 

For the remaining age groups, size msizen is estimated as 



5mc 



8t 



\ mn J \ msizen / 



yrisiz. 



(12) 



Therefore, the size in age groups 2-9 only depends on 
the number of mosquitoes surviving from one age group 
to the next [mn-i) and the size of mosquitoes in the 
younger age group (msizen-i)- 

Model forcing 

To drive a dynamic malaria model it is necessary to have 
boundary conditions that are consistent over time and 
space. Temperature, relative humidity, and rainfall data 
from weather stations are point measures. Hence, they 
might not be representative of larger areas over shorter 



time scales. This is especially true in areas with varying 
topography or where convective rainfall is dominant 
[73-75]. Despite the limitations of rainfall stations, they 
can provide a robust estimate of large-scale events. By 
pooling data from several stations, the error for a sin- 
gle station is reduced and the data can provide a good 
estimate for dry and wet years, for example. Hence, 
weather stations are useful tools for validating climate 
models. 

The problems of point measurements are described 
later, and represent one of the reasons why OMaWa is 
tightly linked to a climate model. As shown in sensitiv- 
ity experiments, the model can also be run with con- 
stant forcing (e.g. temperature) or with data from weather 
stations. 

Where we present results for Africa as a whole, OMaWa 
is driven by data from WRF 3.3.1. This realization (TC50), 
described in part two of this paper, has a tropical channel 
set-up in which set-up, the domain consists of bound- 
aries above and below a certain latitude and no side 
boundaries. The model was run at 50-km resolution from 
January 1, 1989 to January 1, 2009. At the northern (45°N) 
and southern (-45°N) boundaries the model was driven 
by Era Interim. The Kain Frisch cumulus parametriza- 
tion scheme was used [76,77]. This experiment was not 
designed to reproduce observed year-to-year weather 
variability, but to assess the mean mosquito density and 
distribution. The driving experiment is described in the 
section on model validation. 

Climate and weather models 

Currently, our best guess of (future) climate at mul- 
tidecadal time scales comes from general circula- 
tion models (GCMs). These models are designed to 
close the energy budget of the Earth and include an 
interactive representation of the atmosphere, ocean, land, 
and sea ice. A set of scenarios with different emissions 
describes how sensitive the climate is to atmospheric 
constituents (greenhouse gasses) [78]. While climate is 
the average weather over time and space, weather can 
change over minutes, hours, days and seasons. The same 
equations used to predict climate are used to predict 
weather. However, weather forecasts are more depen- 
dent on current observations of the atmosphere. Hence, 
weather predictions are initial value problems, whereas 
climate simulations are rather boundary value problems. 

Both climate and weather models are mostly structured 
on a grid, with coordinates from west to east (x), north to 
south iy) and bottom to top (z). In the grid, one square 
(or polygon) represents the weather within that square. 
While climate models often have a horizontal resolution 
of more than 10000 km^, operational weather models such 
as the European Centre for Medium-Range Weather Fore- 
cast (ECMWF) model are run at approximately 160 krn^. 
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If the state of the atmosphere is observed correctly, higher 
resolution can lead to better local skill in predicting the 
weather. A hybrid between a weather model and a climate 
model is a limited-area model (LAM), which relies on 
initial and boundary conditions from a weather or cli- 
mate model. Given these conditions (weather), the LAM 
can be run at a higher resolution over a limited area, 
which potentially improves the spatial accuracy of the 
coarse model [79]. The WRF model is a widely used 
LAM [41]. 

In tropical regions, most rainfall comes from convec- 
tive clouds. This type of rainfall is generally intense and 
of short duration. The geographical extent of such rain- 
fall episodes may be limited. Therefore, rainfall measure- 
ments in regions where convective rainfall is dominant 
should be handled with care [74,75,80,81], especially when 
extrapolating station data to areas with no data. While 
station data are accurate at a specific point, climate mod- 
els and satellite estimates give a more general descrip- 
tion of the weather within a certain area; Chen and 
Knutson reviewed how models compare to observations 
at varying scales [82]. Since future climate is projected 
using climate models and considering the limitations 
of weather stations, construction of a mosquito/malaria 
model around a LAM is a good choice. The LAM will 
have higher resolution than most climate models, with 
higher-resolution orography, coastlines, and land use, but 
will still give a general description of the weather within a 
certain area. 

Parametrization schemes in the aquatic stages 

To relate a variable such as mortality to the physical envi- 
ronment, we need simplified equations that describe this 
relationship. An equation in which temperature influ- 
ences mortality only states that there is a relationship 
between the two, but does not explain why temperature 
modifies mortality. In this paper we use parametriza- 
tion schemes to represent the influence of the environ- 
ment on mosquitoes. This section describes the aquatic 
parametrization schemes used, excluding water availabil- 
ity, which is discussed later. 

The aquatic stages comprise eggs, four instar stages, 
and pupae. The number of eggs in a location at any time 
is controlled by the number of potential new eggs laid 
(6), available water (/C), natural and induced mortality 
{Pn/i,e/l/p) and movement from the E to the Li com- 
partment. In addition, 20% instant mortality is introduced 
when rainfall exceeds the infiltration rate. This is in line 
with observations by Paaijmans et al. [83]. The number 
of new eggs is simplified to a function of the number of 
gravid mosquitoes in each age group and their size (mea- 
sured as wing length) based on observations [55,84-86]. 
The critical size is set to a wing length of 2.6 mm, which 
is less than that observed by Lyimo and Takken [85] but 



greater than observations by Yaro et al. [87]. Maximum 
wing length is set to 3.3 mm for An. gambiae s,s, [88,89] 
and 3.7 mm for An. arabiensis [70]. The relationship 
between the number of eggs {e) and wing length 
{y^sizen) is then approximated according to the linear 
relationship 

I (-433.3 -h 166.7 • msizen) ' ^« if^s/^e„ > 2.6 mm 1 
^ I 0 otherwise J' 

(13) 

where is the number of mosquitoes in age group n. 
Note that this limits the number of eggs laid by a sin- 
gle mosquito per gonotrophic cycle to approximately 184, 
which is somewhat less than the number observed by Yaro 
et al. [87], but in line with that reported by Howard et al. 
[90]. 

Estimation of water temperature 

Using the 0-10-cm soil temperature {Tsou) from the 
NOAH land surface model [91-94] to approximate the 
mean water temperature (Twater) in larval habitats, we 
assume that evaporative cooling and heat fluxes at the 
water boundaries are negligible. Hence, the water temper- 
ature is equal to the top soil temperature. Paaijmans et 
al. showed that the 5-cm soil temperature represents the 
water temperature in small ponds reasonably well [95]. 
Therefore, the model will have limited validity in areas 
where larger puddles are the main breeding sites. There is 
also a chance that diurnal fluctuations will be slightly over- 
or underestimated. When a grid cell covers several krn^, 
this effect should be negligible, although we do not have 
data to support this. We hope to improve the prediction of 
water temperature in the future, either by modelling this 
explicitly or using a parametrized version based on data 
from Huang et al. [96]. 

Parametrization of mortality 

We used two approaches to calculate mortality in the 
aquatic stages. In the simpler approach, we assume that 
mortality and development time in the aquatic stages are 
independent of the species. We also assume that the rela- 
tionship between the mortality rate and temperature is the 
same for eggs, instars and pupae. In this method we do not 
consider competition effects as described by Paaijmans 
et al. [8]. This type of parametrization is suitable when 
the model is used for one species only (e.g. if the model 
represents an area where only one of the two species is 
present). 

Species-independent mortality (BLL) 

Data provided by Bayoh and Lindsay [97] were used to 
describe the mortality rate according to Eq. (14) {p < 0.01, 
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= 0.81). We call this the BLL method. Mortality rate 
data are plotted in Figure 2b. 

Pn,L iT^ater) = ( + e'^^'^-^'^^-'A 

\ ^ water I (14) 

where ^A^,L (Tvi/^ter) = PN,E{T^ater) = P N,p(T water) 

the aquatic mortality rate per day and Ty^;ater is the water 
temperature (°C). The constants are given in Table 2. 

Species-dependent mortality (KBLL) 

Kirby et al. reported that the mortality rate of An. gambiae 
s.s. and An, arabiensis is modulated by the presence of 
each other in the temperature range 25 — 35°C [72]. To 
account for this we developed two mortality models, one 
for An. gambiae s.s. and one for An. arabiensis. We call 
this parametrization scheme KBLL. The mortality rates 
are based on data from Bayoh and Lindsay [97] and from 
Kirby et al. [72]. Although Holstein also reported larval 
mortality for [An. gambiae s.s.) when exposed to extreme 
low and high temperatures [98], we did not include these 
data when estimating the mortality curves. However, the 
data are plotted in Figure 2 for comparison. Accord- 
ing to our curves, the An. arabiensis mortality rate will 
increase in the range 25 — 35°C as the relative presence 
of An. gambiae s.s. increases. Conversely, the mortality 
rate of An. gambiae s.s. will decrease as the proportion of 



Table 2 Constants for equation 1 4 and 33 



Constant 


Value 


Equation 


k^ 


700000 


14 


ki 


8.4 


14 


k3 


.126 


14 


kA 


10.8 


14 


ks 


150 


14 


ke 


-.08 


14 


kv 


.1 


14 


k8 


-.61 


14 


k9 


33 


14 


C] 


0.1675256 


33 


C2 


0.0121402 


33 


C3 


0.1686 


33 


C4 


1.991 


33 


C5 


1.881 


33 


ce 


4.641 589e26 


33 


c? 


250 


33 


C8 


23 


33 


Cg 


12 


33 


ClO 


100 


33 




3 


33 
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An, arabiensis increases. The mortality rate ^jqi is given 
by 



Farab = mm I — ^ , 1 I 

\ Z^n=l ^n,gamb / 



(15) 



I 0.002404075 • T^^^^^ - 0.1127944 • Twater + 1.337783 

Pn,l (Twater) ' (0.4 + 0.6 • (1 + sin (-10.9956 + 0.3142 • T^ater)))^'''''' if 25 < Twater < 35 



and 



Fgamb = min 



^n=l ^n,gamb \ 
^^=1 Ln,arab / 



(17) 



,arab twater) 



0.0006556736 • T^^^^^ - 0.02980226 • Twater + 0.3587285 
^N,L iJwater) ' ((2 + COS (-18.8496 + 0.6283 . r^«^.r))^-^^°^°°^)^''"' if 25 < Twater < 35 (18) 

PN,L,gamb (Twater) Twater < 21.91209. 



Fgamb and Farab ^re the ratio of An. gambiae s.s. to An, 
arabiensis larvae and An, arabiensis to An, gambiae s,s, 
larvae, respectively. At each time step, Lsize is estimated as 
a function of Bi and K, As the density increases, there will 
be more competition and hence less food for each larva, 
which leads to smaller larvae. 

Parametrization of the development rate 

The rate of development between the different aquatic 
stages follows the corrected version of Bayoh and Lindsay 
[97]. Since these data are only valid for An, gambiae 5.5., 
we made a small modification to prolong the develop- 
ment times for An, arabiensis. Data from Kirby et al. 
[72] and Paaijmans et al. [8] suggest that time for devel- 
opment from a larva to an adult is approximately 5.5% 
longer for An, arabiensis than for An, gambiae s,s. Hence, 
we increased the development time for An, arabiensis 
by 5.5%. The reason for this longer development time is 
that An, arabiensis takes longer to develop a larger body. 
Curves of the development rate are shown in Figure 3. 

The two previous studies also suggest that the devel- 
opment rate [8] and mortality [72] of the two species 
are modulated by the presence of each other, so we take 
account of this in out model. The development time for 
An, arabiensis is prolonged in the presence of An, gam- 
biae 5.5., while the time is shortened for An, gambiae s,s, as 
the relative proportion of An, arabiensis increases. Using 
data from Paaijmans et al. [8], the development rate r is 
modified according to 



farab = min 100 ■ 



^^n=l ^n,arab 



^^=1 ^n,gamb H" ^2n=l ^n,arab 



-,75 



T^gamb = T^gamb ' (1 - farab ' 0.0008421) ^ 

for An, gambiae s,s, and 



(20) 



fgamb = min 100 



^n=l ^n,gamb 



^^=1 ^n,gamb ~l~ ^2n=l ^n,arab 



-,75 



(19) 



(21) 



^arab — ^arab ' 

(1 • 0.002138) (22) 
for An, arabiensis, farab fgamb is the fraction of An, 
arabiensis and A«. gambiae s,s,, respectively. 

Parametrization of breeding sites 

The formation of puddles can be described as a balance 
of runoff, infiltration, evaporation, and rainfall entering 
the puddle. The formulation of an idealized puddle can be 
found in Additional file 2. 

Modelling of every single breeding site requires high 
enough resolution to resolve the puddle. In practice this is 
not possible and the problem has to be simplified. 

Mushinzimana et al. described typical breeding sites in 
a Kenyan highland area [99]. Most of the puddles were 
located at less than 100 m from rivers, which means 
we can assume that semi-permanent puddles will mostly 
form in the proximity of rivers and lakes. They also found 
that the number of breeding sites was close to threefold 
higher in the rainy season compared to the dry season, and 
grouped breeding sites by surface area. 

If we assume that breeding mainly occurs in the vicin- 
ity of potential rivers and lakes, the availability of breeding 
sites can be expressed as a function of potential river 
length and soil saturation. At high resolution this might 
not always be true [6], but since the model is designed 
to be applied to coarser grids, we believe the assumption 
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Figure 3 Water temperature according to development time in days from first instar to adult. Left panel: ratio of An. gambioes.s. to An. 
arabiensis. When greater numbers of An. gambiae s.s. are present, An. arabiensis develop more slowly. Right panel: ratio of An. arabiensis to An. 
gambiae S.S.. When greater numbers of An. arabiensis are present, /\n. gambiae s.s. develop more quickly. 



is as reasonable as or more reasonable than the common 
assumption that puddle formation is only dependent on 
rainfall [29]. The newest version of the NOAH land sur- 
face model in WRF 3.4 also includes groundwater and 
dynamic vegetation, and future versions might change the 
way in which puddles are parametrized. In OMaWa we 
introduce a simple parametrization scheme to represent 
breeding sites. 

The Hydrological Data and Maps based on SHuttle 
Elevation Derivatives at Multiple Scales (HydroSHEDS) 
15s river data set from the US Geological Survey (USGS) 
[100] was used to derive the total potential river length 
within a grid cell. Since the algorithm used to develop this 
data set describes where water would collect if it were 
available within the catchment, it also represents a general 
description of the potential for water aggregation within 
an area. However, the validity might decrease on moving 
to finer scales [6]. 

Here we divide rivers into three different classes: peren- 
nial, intermittent and ephemeral streams. For each class, 
potential river length {Rp, km) within a grid is defined as 



27TER 
360 



• cos (p, 



(23) 



where S is the equally spaced river data-set resolution in 
degrees, where Mon = Alat, ER is the radius of the Earth 
(6371.22 km) and (p is latitude in radians. 

In a simplified model we estimate puddle volume as 
a function of river length and relative soil moisture. 
Although this is a very crude estimate, we compared this 



simple model with data from Mushinzimana et al. [99] and 
derived a simple expression for the carrying capacity in a 
grid cell: 



K = . . SMr 

KWl river 



(24) 



where is the maximum larval biomass per km of 

Km river 

river (2400 mg, estimated from data collected by Munga 
et al. [101]) and SMr is the relative soil moisture content 
(fraction). 

In the current implementation we do not distinguish 
between fast- and slow-flowing rivers. It should be noted 
that this way of approximating breeding sites has limited 
validity in areas with irrigation or around rivers where 
breeding sites could form as rivers recede [66,67,102]. 
Some special cases, such as along the River Nile in Sudan, 
where breeding sites form as a result of rainfall hundreds 
of kilometers away, will not be captured at all [103]. 

Parametrization of the gonotrophic cycle 

The gonotrophic cycle depends on temperature and 
is important for the vectorial capacity of mosquitoes. 
Lardeux et al. studied the gonotrophic cycle for An. pseu- 
dopunctipennis [104]. We combine their data with other 
published studies on anophelines to estimate the length 
of the gonotropic cycle. There are few studies on An. 
gambiae s.L, and hence we have to assume that other 
anophelines share the same physiology and strategy with 
respect to the gonotropic cycle. Ruiz et al. showed there 
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are some differences [23], but until further evidence of the 
reproductive strategies of different members of Anophe- 
les genera becomes available, we will not consider this 
effect. Studies used to develop the formula include those 
by Guillermo et al ([105], An, alhimanus), Afrane et al 
([106], An, gamhiae s,l,), and Maharaj ([107], An, arabi- 
ensis). We also include the formula given by Hoshen and 
Morse [108]. Their model is based on degree days and 
is included according to Eq. (26). The gonotropic rate 
{day~^) and data used to develop the formula are shown 
in Figure 4. 

Fgonot = min (max ( - ^ + ^ • ^^/r, 0 ) , .5 ) (25) 



G{T) = U + — — — ) 

' Fgonot + (l.71 + 544347.6 • tJ-^^Y^^ • (l - F^onot) > 

(26) 

where Tat is the air temperature (°C), is degree days, 
and Tc is the critical temperature from Hoshen and Morse 
[108], withDj = 37, and Tc = 7.7. 



Parametrization of the age-dependent mortality of adult 
mosquitoes 

The mortality of adult anophelines differs according to 
age and species [7,107,109]. This has often been over- 
looked in mosquito models [23,110]. To show how this 
assumption can influence the stability of mosquito pop- 
ulations and malaria transmission, we use the mortality 
model of Martens [110] as a reference. We also plot Eq. 
7 from Ermert [29] in Figure 5 to highlight the differ- 
ences between this model and established models. For 
convenience, we repeated Martens equation, as follows: 



pM,m(T) = l-e -4.4+i.3i.r-.o3.r2 ^ (27) 

Our new survival curves are based on unpublished data 
from Bayoh and Lindsay [47]. The validity ranges from 
5 to 40°C by 5°C and 40 - 100% by 20% relative humid- 
ity. We name the scheme BLLad (Bayoh-Lindsay-Lunde 
adult mortality). The data set and the curves are valid 
for An, gambiae s,s. The lowest agreement between the 
model and the data is at 40% relative humidity and 40°C. 
While the data suggest that all An, gambiae s,s, would 
be dead after approximately 2 days, the survival curve 
would result in no mosquitoes after approximately 4 days 
at 40% relative humidity and 40°C. To correct for this 
error, we include data from Kirby and Lindsay [111], who 
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Figure 4 Inverse of the duration of the gonotropic cycle according to the mean daily temperature (in °C). The solid black line shows Eq 26 
and the dashed line shows the formula given by Hoshen and Morse [1 08]. 
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Figure 5 Proportion of An. gambiaes.s. surviving at 60% relative humidity and mean temperature of 0, 10, 20, 30, and 40°C (selected for 
clarity) according to time (in days). Dashed vertical lines indicate the different age groups in the model. The Survival curve panel shows Eq. 31, 
while the Numerical solution panel shows survival in the model when the age groups are split into nine classes. For reference we also show survival 
according to the Marten equation (27) and Eq. 7 from Ermert et al. (Bayoh scheme) [29]. The mean absolute error for all combinations of 
temperature and relative humidity was 73 for our model, 1 71 for the Marten model, and 1 29 for the Bayoh scheme. 



described the responses of An. gambiae s.s. and An. ara- 
biensis to high temperatures. By assuming that maximum 
survival is 480 min for An. gambiae s.s. and 1440 min 
for An. arabiensis at temperatures greater than 40°C, we 
can set the mortality rate to 3day~^ and lday~^y inde- 
pendent of age group. However, there are uncertainties at 
relative humidity below 40%. The lack of studies in this 
range is a limitation of this survival model, and could 
make the model less accurate for An. gambiae s.l. in 
some regions. The basic principle of these survival curves 
is that mortality will be low in the first few days after 
emergence. In addition, mosquitoes that survive up to a 
certain age have a higher survival probability (depend- 
ing on Tair and relative humidity). In Figure 5, survival 
at 60% relative humidity and 0, 10, 20, 30, and 40°C 
is plotted. 

Size affects the survival of adult mosquitoes [7,51,54,55] 
([56], Aedes aegypti). If we assume that the major differ- 
ences in mortality between An. gambiae s.s. and An. ara- 
biensis can be attributed to mosquito size, we can modify 
a as a linear function of mosquito size. Here we subjec- 
tively choose reasonable constants for h (msize)- T^air may 
be completely or partly replaced by indoor temperature 
{Tmdoor> described later), depending on the proportion of 
mosquitoes indoors. In experiments covering the African 
domain, we assumed that 80% of An. gambiae s.s. and 20% 
oiAn. arabiensis are located indoors. 

g {nisizen) = 2.1731 - 0.3846 • msize (28) 



f(RH) = 6.48007 + 0.69570 • (1 - e'^-^^'^^) (29) 



01 =g (msize) 

^lo+(l+^)^'/'^((l+^)^(l+^).2-/(i?//)) 



X e 



(30) 



mM,m(oi,^,a) = ^ 



i=0 



I (a ' a)^i=o 



1 



J-a-a) 



(31) 



where f = 6, ^ is a function of mosquito size, and RH is 
relative humidity. The mortality rate for each age interval 
can then be approximated as 



log 



At 



ifr < 40 
otherwise 



(32) 



If we assume that differences in adult mortality for An. 
gambiae s.s. and An. arabiensis can be explained by dif- 
ferences in body size, these BLLad curves can be used for 
both species. We explore this mortality model in [64]. 

AL adult mortality 

A similar approach can be used for An. arabiensis. Using 
survival curves reported by Afrane et al. ([112], Figure 
two) (copyedited with g3data [113]), we can estimate mor- 
tality based on the daily maximum temperature. Because 
of the few data points, this approach is much more 
uncertain and should be considered experimental. The 
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advantage of this mortality model is that the data are not 
estimated from a laboratory setting. The maximum tem- 
perature reflects some aspects, such as radiation, albedo, 
and humidity, of the environment in which mosquitoes 
live. In some of the results presented in part two, we use 
this model for adult survival. 



Oi = Ci - C2 
I 



° mod V mod 



- Tmod ■ C5 - Cs) - C7 • 



r 

C6 



+ 18 



11.10 



(33) 



(34) 



Constants ci,...,ii are listed in Table 2. By setting ^ = 2 
we can simplify the survival curve for An, arabiensis to 



(35) 



The corresponding curve is shown in Figure 6. 



Parametrization of air temperature 

Paaijmans et al. discussed the importance of using indoor 
rather than outdoor temperature, to describe the environ- 
ment for mosquitoes and parasites [114]. They included 
two studies that showed the relationship between indoor 
and outdoor temperature in Kenya [115] and Tanzania 
[116]. Here we add two additional studies, one from Kenya 
[48] and one describing the temperature in traditional 
and low-cost modern housing in the Eastern Cape, South 
Africa [117]. The data used to parametrize equation 36 
came from; 1, Afrane et al. [48]; 2, Makaka and Meyer 
[117]; and 3, Paaijmans et al. [114-116] {R^ = 0.89). It 
is clear that temperatures inside a house are more stable 
than outdoor temperatures. House type greatly influences 
daily temperature fluctuations [117,118], and the model 
used here might not be valid for all house types. While 
some studies have assumed that houses are always hotter 
than the surroundings [119], we approximate the indoor 
temperature as 



Tindoor = 10.33 + 0.58 . Ta, 



(36) 



Since the data are based on maximum and minimum 
temperatures, the timing of the indoor temperature might 
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Figure 6 Proportion of An. arabiensis surviving at daily maximum temperatures. Estimated from Afrane et al. [1 1 2] (blue line). Dashed vertical 
lines indicate the different age groups in the model (grey lines). 
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be offset by a couple of hours. This is evident in a study 
by Makaka and Meyer [117], who delayed the maximum 
indoor temperature by a couple of hours compared to the 
environmental temperature. At present we do not account 
for this delay, since the diurnal temperature ranges will be 
correct even if we do not. The data and regression line are 
shown in Figure 7. Further studies on indoor compared to 
outdoor temperatures are needed to make this correction 
more accurate. 

Hence, Tat can be partly or fully replaced by Ti^door> 
depending on the proportion of mosquitoes indoors. 

It should be noted that we still do not include tem- 
peratures in resting places described by Holstein, such as 
holes in rocks and cracks in soil, covered pigsties, rabbit 
hutches, hen coops and dry wells [98], and by de Meillon 
([120], under stones). 

Approximation of mosquito movement 

The role of diffusion and advection in vector borne dis- 
eases have been explored in several papers [102,121-127]. 
Considering the gradual invasion of Brazil in the 1930s 
by An. arabiensis [60] it can be argued that move- 
ment of mosquitoes is important over decades. Here we 
include the active and passive transport of mosquitoes 
as fluxes across grid boundaries. Passive transport is 



movement of mosquitoes caused by wind, while active 
transport is movement due to flying. On shorter time 
scales the role of such movement will be limited. 
However, on long time scales it is necessary to allow 
mosquitoes to travel to allow them to establish in new 
locations. 

Transport of mosquitoes is defined by fluxes {s~^) 
at the grid boundaries. In the model we allow fluxes 
from the eight neighboring grid points. A special case 
is implemented when a neighbouring cell is water. In 
this case, fluxes to water are reduced to 0.1% of the 
original flux to avoid large losses of mosquitoes along 
the coastline. Given strong winds from land to the 
ocean, such an assumption could lead to accumulation 
of mosquitoes along the coast. Conversely, allowing free 
movement to the ocean could lead to undesired loss of 
mosquitoes. 

Since the movement of mosquitoes has a high com- 
putational cost, the spatial fluxes do not change the 
size calculations. This will introduce some minor errors 
when the movement of mosquitoes is low compared 
to their density, with larger errors if many mosquitoes 
are moved relative to their density. When a cell 
free of mosquitoes is colonized, the size is set to 
3.05 mm. 



40- 




10- 



10 20 30 40 

Outdoor temperature 

Figure 7 Relationship between outdoor and indoor temperatures. Numbers denote the study from which data were taken: 1 , Afrane et al. [48]; 
2, Makaka and Meyer [1 1 7]; and 3, Paaijmans et al. [1 1 4-1 1 6]. The blue area represents the 95% confidence interval, and the black line shows Eq. 36. 
/?2 = 0.89. 
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The possible flight range of anophelines varies with food 
availability [128]. We do not include vegetation types in 
the model and hence it is hard to justify differences in 
flight performance based on, for example, land use. The 
dispersion coefficient describes how far mosquitoes can 
move in a day. We assume that the dispersion coeffi- 
cient D is constant, independent of geographical location. 
For An. gambiae s.s. and An. arabiensis, real flight per- 
formance outside the laboratory of only a few hundred 
meters per day (approx. 300-700 m) has been reported 
[102,129,130]. In this experiment we subjectively chose 
D = 30mday~^ independent of age group. Anopheli- 
nae also travel with humans [131], which adds to the 
transport equation and makes the dispersion coefficient 
uncertain. Gillies noted that wind direction mostly has a 
minor effect on dispersal [129], while de Meillon [132] and 
Adams [133] reported distances of 2-4.5 miles (3-7 km) 
in the direction of the prevailing wind. Thus, it can- 
not be ruled out that wind plays a role on longer time 
scales. Hence, we express movement caused by wind as 
a function of 10-m zonal (u) and meridional (v) wind 
components {ms~^). This can be understood by consid- 
ering the following example. For a constant r/-wind of 
10ms~^ and v-wind set to 0, mosquitoes will be moved 
a distance related to a scale factor Sf, which is equal to 
the distance travelled at 20ms~^ to the east. For exam- 
ple, with Sf = 750mday~^y the eastward distance traveled 

will be Sf ' 2Qms-^ ^"^^ m in 1 day, but since each mosquito 
is not modelled individually, it would be more natural 
to describe this as a fraction moving a certain distance. 
Different wind directions and speeds will result in other 
distances/fractions and directions. D and Sf are unknown 
tunable constants. 

Since the species considered here are most active at 
night [22], movement will be suppressed between 06:00 
and 18:00 h (local time) and amplified at night according 
to 



(cos(ir.^) + i)' 



1.506925 



(37) 



where LT is local time, f^^ k ^ 1 and 



LT= 



UTCtime + 



15 

longitude 



15 

longitude 
15 



+24 if^rQ/^,+ 



15 

longitude 
15 



<=0 



Transport of mosquitoes and mosquito sizes inside and 
outside a grid are defined by 



^ = E E ^^■J^n,>J- (38) 

i=-ij=-i 

More specifically, during a time A^, movement can be 
calculated as follows. On a day with no wind, transport is 
equal in all directions, D = 30mday~^, and the flux at a 
boundary is defined as 



D 



Adij . 24 • 60 • 60 



,^ ={-l,l},7 ={-1,1} 



(39) 



and transport rjij^ is then equal to 



otherwise. 



In the presence of wind, we obtain additional trans- 
port as a function of zonal and meridional wind 
components. 

Mortality related to feeding 

One factor that is often overlooked in malaria (mosquito) 
models is survival related to food availability {P(B)). 
Ye-Ebiyo et al. reported that maize pollen availability has 
a positive effect on larval (and hence mosquito) fitness 
[70,134]. Creating maps of plant types is beyond the scope 
of this study, and hence we chose not to account for 
mortality related to crops. However, we performed initial 
tests in which we included GlobCover Land Cover version 
V2.2 (European Space Agency [135]) to give a rough esti- 
mate of regions where increased fitness could be expected. 
The other source of food for female anophelines is blood. 
Compared to a starved mosquito, a mosquito that has 
had access to blood on days 1-3 has a theoretical flight 
distance that is increased by a factor of 6-7 [128]. There- 
fore, it is plausible that the higher (lower) the probability 
of finding a blood meal {P(B))y the higher (lower) is sur- 
vival in the early life stages of adult mosquitoes. Bouma 
and Rowland reported higher parasite prevalence among 
children of families who kept cattle compared to those 
who did not [136], which can indicate either higher sur- 
vival (older mosquitoes) or simply that some anophelines 
are attracted to cattle. If we assume that a newly emerged 
mosquito has a flight range of Fr^ = 0.5km^day~^, the 
daily probability of finding a blood meal can be calculated 



P(B) = 



HBI • Phuman ' + (1 - HBI) • Pbovine ' Pryn if P{B) < 1 ] 

1 Otherwise I 



(41) 
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where Phuman ^rid Pbovine is the probabiUty of finding a 
human and bovine source, respectively. Phumans is defined 
as the human population density per knP' multiplied by 0.1 
(since a smaller area on a human is accessible) and Pbovine 
is defined as the bovine density per knP'y each with a user- 
defined threshold at which the density is so low that P{B) 
is virtually zero. Since P{B) is a conceptual parameter, it 
can be tuned. 

Since blood meals, besides sugar meals, are important 
for the mobility [128] and survival of female anophe- 
lines [137], the success of a species is likely to be linked 
to the presence of the preferred host. The dominant 
blood source for An, arabiensis is bovine and human 
blood, while it is human blood for An, gambiae s,s, [138]. 
In reality there are strong indications that the human 
blood index is a dynamic quantity rather than a constant 
[139-142]. In the current implementation, HBI is a static 
number and hence there are probably errors related to this 
term. To find the probability of feeding on humans at each 
time step, we combine two data sets. Between 2000 and 
2010 we use population densities from the Gridded Popu- 
lation of the World (GPW) [42], and for before 2000 and 
after 2010 we use growth rates from the Population Divi- 
sion of the Department of Economic, and Social Affairs 
of the United Nations Secretariat [143]. Since there are 
no projections of cattle densities, this quantity is time- 
invariant and based on Food and Agriculture Organiza- 
tion (FAO) 2005 estimates [44]. We are currently working 
to include time-varying cattle densities. 

In the model, mortality caused by food limitations is 
a function of how many humans or cattle are available 
per mosquito and the human blood index. We assume 
that HBI is time- and space-invariant, and only depends 
on the species. For simplicity we chose available humans 
to be humans who are not sleeping under a bed net. In 
the simulations presented here, we set bed net usage to 
zero, and hence the results represent mosquito distribu- 
tion without interventions. Bayoh et al. hypothesized that 
the survival of the different species is related to the avail- 
ability of the preferred host [9]. The daily mortality rate 
caused by limited human blood is expressed as 

/ -^Th \ 
Ph,m = max[l- ™^ ,0 . (42) 

The functional form of of equation 42 can be seen in 
Additional file 3. 

Figure 8 shows the probability of finding a blood meal 
for the sibling species on January 1, 1999. 

Results and discussion 

Sensitivity experiments 

Sensitivity experiments are useful in understanding which 
parameters are important for the success oiAn, arabiensis 



and An, gambiae s,s, and which are important for malaria 
transmission. Classical sensitivity analysis investigates the 
robustness of a study when parameters are estimated from 
statistical modelling. Our model uses parametrization 
schemes to represent the influence of the environment 
on the two species. We show how the model responds 
to changing temperature, humidity, mosquito size, disper- 
sion and the probability of finding blood. This approach 
does not allow us to directly measure the robustness of 
each parametrization scheme, but gives us an insight into 
which external factors influence the model and where 
it is of importance to have improved parametrization 
schemes. We use the term sensitivity experiments for this 
analysis. 

Settings 

To demonstrate some of the capabilities of the model, we 
set up a series of experiments. Some aspects are best visu- 
alized as a one-dimensional model (time and age), while 
other features are shown using a spatial domain (time, 
age, and space). For the one-dimensional experiments, the 
water temperature is set to the air temperature, except for 
temperature greater than 33°C, for which we set temper- 
ature to 33°C. This modification is required since pupae 
and fourth instar larvae will not develop below 18°C or 
above 34°C [144]. The results are therefore less robust 
when temperature is greater than 33°C. Unless other- 
wise stated, we use size-dependent mortality, correction 
for indoor temperature, the KBLL method to estimate 
mortality in the aquatic stages, correction for the develop- 
ment rate in the aquatic stages depending on the ratio of 
each species, and movement of mosquitoes (in the spatial 
cases). 

Sensitivity to temperature, relative humidity and 
mosquito size (TempHumSize) 

The age-dependent mortality is influenced by tempera- 
ture, relative humidity and mosquito size [Eq. (32)]. This 
experiment explores how the dynamics of malaria is sen- 
sitive to temperature, relative humidity and mosquito size 
(measured as mm). We assume that no births occur to 
isolate the effect of the transmission process, and con- 
sequently constant mosquito body size in the course of 
integration, but include mortality and the biting rate. 
In this experiment we assume that only one species is 
present (since the main competition occurs in the aquatic 
stages). This sensitivity test is designed to observe how 
the proportion of mosquitoes becomes infected as a func- 
tion of temperature, relative humidity and mosquito size, 
given that we start with 1000 newly emerged mosquitoes, 
with mi = 1000 and m2-9 = 0 as the initial condi- 
tions. In this experiment, 1% of the human population 
is infectious for Plasmodium falciparum. Mosquitoes are 
infected with an efficiency of 100%, meaning that biting 
an infectious human results in gametocyte transmission to 
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Figure 8 Probability of finding a blood meal for >An. orobiensis [HBI = 0.4) and An. gambiae s.s. [HBI = 0.95) with zero bed net coverage. 



the mosquito. In practice, this would be the same as say- 
ing that 10% of humans were infectious and gametocyte 
transmission had an efficiency of 10%. We also neglect 
the effect of heterogeneous biting. This is the only exper- 
iment in which we model the proportion of infectious 
mosquitoes explicitly. The modified equations describing 
the transmission process are described in [64]. 

The rate of sporozoite development within mosquitoes 
is expressed as [5] 



pf =[a^ 



(43) 



where a = 9.5907, b = 0.0051029, c = 0.7349, and d = 
17.0325. This expression was derived from the figure in 
MacDonald page 119 [5] using g3data [113], and fitted 
using non-linear least-squares [145]. 

The gonotrophic cycle and biting rate are defined in Eq. 
(26). 

The integrations are repeated with different combi- 
nations of temperature and relative humidity. This is 
a simple representation of gametocyte transmission to 
mosquitoes and is an idealized approach for exploring 
the proportion of mosquitoes (of the original 1000) that 
would become infected under different temperature, RH 
and mosquito size. Figure 9 shows how the percent- 
age of infectious mosquitoes changes with temperature, 
RH and mosquito size. Lyimo and Koella reported that 
the largest mosquitoes were less likely to have sporo- 
zoites, but had more oocysts than smaller mosquitoes 



[54]. They attributed this to increased mortality in the 
presence of many oocysts, an effect that is not included 
in our model. Figure 9 shows that the potential per- 
centage of infected mosquitoes is sensitive to all three 
parameters in the model. Although higher survival has 
been attributed to body size in dry [7,51] and semi-arid 
environments [55], the advantage or disadvantage of a 
larger body has been poorly described in saturated envi- 
ronments. Therefore, the sensitivity to body size at 80% 
RH should be interpreted with care. According to the 
model, temperature is not the only factor that governs the 
transmission of malaria (in areas with no interventions); 
humidity and how mosquitoes adapt to dehydration stress 
are also important factors. The most efficient transmis- 
sion, expressed as the integral, with respect to days, occurs 
at 25°C at 40% and 80% RH, and at 24.5°C at 10% RH, 
independent of mosquito size. 

These results should be viewed in light of recent find- 
ings by Paaijmans et al. that optimal transmission occurs 
at lower temperatures [4]. 

Sensitivity to temperature and carrying capacity 
(TempCar) 

The aim of this sensitivity test was to investigate how 
carrying capacity and temperature determine the relative 
proportion of An. arabiensis and An. gambiae s.s. We 
set the relative humidity to 80% and the probability of 
getting a blood meal to one. We assumed that the soil was 
saturated and we varied the temperature between 16 and 
38°C (with corrections over 33°C for water temperature) 
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Figure 9 Percentage of 1000 mosquitoes that are infectious after x days. The y-axis represents temperature in degrees centigrade. The model 
is integrated at two mosquito sizes (2.8 and 3.2 mm for wing length) and three relative humidity values. 



and the carrying capacity between 0.0625 and 
125 mgkm~^. 

Carrying capacity in the aquatic stages influences larval 
growth and adult survival. While An. arabiensis invests 
more time in growth than An. gamhiae s.s., the former 
develops a larger body, and consequently has the poten- 
tial to oviposit more eggs than the latter. If the two species 
experience the same mortality rate in the aquatic stages, 
more An. gambiae s.s. will emerge, but over time An. ara- 
biensis can face this challenge by outnumbering the eggs 
of An. gambiae s.s. in the habitat. Thus, we are interested 
in testing how the carrying capacity in the aquatic stages 
alters the relative proportion of each of the adult species. 
In this model we only consider the competition between 
these two species, and hence neglect other competing 
species [146]. 

As observed in Figure 10, An. gambiae s.s. dominates 
between 27 and 30°C. This is the eff'ect of the develop- 
ment rate modifications described by Kirby et al. [72] 
and Paaijmans et al. [8] (Figure 2 and "Species-dependent 
mortality (KBLL)"). Interestingly, the dominance of An. 
arabiensis is most pronounced in the drier simulations, 
meaning that high competition, compared to adult sur- 
vival, is favourable for this species. This can be attributed 
to the strategy of larger body size and higher egg produc- 
tion. Lehmann et al. found that An. arabiensis dominated 
during the dry season, while An. gambiae s.s. dominated in 
the rainy season [57]. The advantage oi An. arabiensis in 
crowded breeding places might be one factor contributing 



to the shift in species composition as the surface area of 
puddles starts to shrink. 

Sensitivity to temperature and the probability of finding 
blood (pBloodlD) 

This experiment shows how the model responds to 
changes in the probability of finding a blood meal, which 
influences the rate at which mosquitoes can oviposit and 
increases energy consumption if hosts are hard to locate. 
If, for example, cattle are easier to find compared to 
humans. An. arabiensis will potentially use less energy 
per batch of eggs and will also be able to utilize breeding 
sites at a higher rate than An. gambiae s.s. It is also pos- 
sible that An, arabiensis uses cattle for navigation [147]. 
Over time, such differences might lead to dominance by 
one species. In this experiment, we varied the probabil- 
ity of finding blood, P{B), for An. arabiensis from zero to 
one, as well as varying the temperature as described for 
TempCar. 

We set the probability of finding blood to one for An. 
gambiae s.s., independent of the probability of An. ara- 
biensis finding a blood meal. This is a purely theoretical 
experiment designed to demonstrate a concept. The prob- 
ability of finding blood is varied between zero and one for 
An. arabiensis. The scenario in which P(E) = 1 for An. 
gambiae s.s. and zero for An. arabiensis is not a realistic 
scenario, but the difference in P(B) is grounded in differ- 
ences in their feeding behaviour, whereby An. arabiensis 
can utilize cattle more efficiently than gambiae s.s., for 
example. 
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Figure 10 Fraction of ^n. arabiensis as a function of air temperature and carrying capacity. The water temperature is set to the same value as 
the air temperature unless the temperature is greater than 33°C (at which most pupae would not develop into adults [1 44]). In this case the water 
temperature is set to 33°C, but the productivity will remain low. The fraction of An. gombiaes.s. is one minus the fraction for An. arabiensis. 



Figure 11 shows the relative fraction of An. arabiensis. 
In addition to the pattern observed in Figure 10, it is also 
evident that if P{B) is low for An, arabiensis, An, gambiae 
s,s, dominates. P{B) can be interpreted as a parameter that 
describes both the probability of finding blood for repro- 
duction and survival, and the energy spent in the search 
for a blood meal. For example, easy access to cattle might 
give An, arabiensis an advantage in exploiting breeding 
sites, which could lead to suppression of the number of 
An, gambiae s,s, if increased use of bed nets reduces the 
effective human population or causes higher mortality of 
anthropophilic species. This mechanism might help to 
explain the decline in An, gambiae s,s, observed by Bayoh 
et al. [9]. 

Sensitivity to the probability of finding blood in a spatial 
domain (pBloodlD) 

This experiment is similar to pBloodlD, but this time we 
integrate the model for 5 years over the African domain. 
The experiment consists of two runs, for which the first 
has P{B) similar to Figure 8 and the second has P{B) = 1 



over all land areas for both species. The population den- 
sity is space-invariant at 400 humans /km^ (remember that 
the number of mosquitoes is limited by the number of 
hosts). Thus, the only limitation in this experiment is the 
physical environment (air and water temperatures, rela- 
tive humidity, wind and run-off), which is updated every 
3 h. The initial conditions for the mosquito populations 
were the same for the two runs. 

Even though we have stated that the probability of find- 
ing blood P{B) is an expression of the cost of finding a 
host, it might well be that P{B) also includes a compo- 
nent that describes the environment shaped by cattle and 
humans. Therefore, it should be noted that it is difficult to 
distinguish between the true probability of finding blood 
and the environmental changes caused by the presence of 
humans or cattle. 

Under the scenario of equal probability of finding blood 
for the two species. An, gambiae s,s, loses the competi- 
tion after 5 years (Figures 12 and 13), probably because of 
the greater reproductive potential of An, arabiensis. The 
only strongholds left for this species are DRC, Congo, and 
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Figure 11 Fraction of An. arabiensis as a function of air temperature and probability of finding a blood meal. The fraction of An. gambiaes.s. 
is one minus tine fraction for An. arabiensis. 
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Figure 12 Square root of number of ^n. arabiensis per km^ in the two pBlood2D experiments. In PO we used realistic values of the probability 
of finding blood, P(B), while in PI the probability of finding blood was set to 1, independent of the location. See the text for further details. 
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Gabon. Hence, the strategy of An. arabiensis to develop 
a larger body, produce more eggs, and possibly reduce 
adult mortality at the cost of spending more time in the 
aquatic stages is successful when access to blood is unlim- 
ited. An, arabiensis has extended its distribution as far 
north as the southern tip of Western Sahara. While the 
original set-up of the model (PO) limits the distribution 
of An. gambiae s.l. to approximately 17°N in the Sahel, 
the experiment with P{B) = 1 (PI) has a distribution up 
to 22°N in Mali, Niger, Chad and Sudan. This is in line 
with observations of the northerly limit of An. gambiae 
s.l. [148-150]. The lack of An. gambiae s.l. north of 17° in 
the original set-up (PO) might be a result of the way the 
model is formulated. The population density is calculated 
within a box of approximately 50 km x 50 km. It might 
well be the case that pockets of higher population/cattle 
densities within this box could sustain a mosquito popu- 
lation. This is not resolved in the model. It is also worth 
mentioning the study by De Meillon [151] of the anophe- 
lines of Namibia, which revealed that An. gambiae s.l. is 
present in large parts of the country. The original set-up 
(PO) allows sustainable mosquito populations in Namibia, 
while the density of An. gambiae s.l. in PI is more compa- 
rable to the observations of De Meillon. The problems of 
capturing the distribution of An. gambiae s.l. in Namibia 
may originate from the problems of resolving pockets of 
high host density or changes in cattle density and distri- 
bution at the time of the study compared to the present 
day [43,44,152]. 



It is also worth mentioning that the density of An. 
gambiae s.l. in South Africa is not very sensitive to the 
probability of finding a blood meal. Hence, the distribu- 
tion of An. gambiae s.l. is mainly restricted by climate 
according to the model. 

Figures 12 and 13 show the distribution and density of 
An. arabiensis and An. gambiae s.s. under realistic (PO) 
and space-invariant (PI) P{B) after 6, 12, and 18 months. 
The integration was started on January 1 and the model 
was run for 5 years. 

Mosquito transport (mosqTran) 

The purpose of this experiment was to demonstrate 
how the initial conditions and competition influence the 
distribution of An. gambiae s.s. and An. arabiensis. To 
explore the theoretical dispersion distance and the influ- 
ence of the initial conditions, we set up a simple exper- 
iment. In mosqTran(a) the model was initialized with 
An. arabiensis at -4.494381°E, 14.0154°N (Sahel), and 
An. gambiae s.s. at -4.494381°E, 6.502846°N (Cte dlvore. 
Ivory Coast) on January 1, 1989. The second experi- 
ment, mosqTran(b), had the same setup, but without 
An. arabiensis. 

The purpose of this demonstration was to show the 
importance of mosquito movement and how new areas 
can or cannot be colonized. In a model in which move- 
ment is restricted, the vector range would also be 
restricted by the initial model conditions. For example, 
if only one point was specified for mosquitoes at the 
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beginning of the integration, only the same point would 
have mosquitoes after 100 years. With dynamic move- 
ment the mosquitoes could colonize new areas if the 
environmental conditions, or the probability of finding 
blood, change over time. 

Figure 14 shows the relative difference in An. gambiae 
s,s, distribution in the two experiments. It is evident that 
in the presence of An, arabiensis, An. gambiae s.s. fails to 
colonize large parts of Mali and Burkina Faso. It can be 
argued that this is not a result of the initial conditions, 
but of competition. Additional file 4 illustrates why this 
is indeed a result of the initial conditions, although the 
initial conditions would not play a role in the absence of 
competition. 

Figure 15 shows the number of months required to 
reach a density of 20 mosquitoes / kn? . It is interesting to 
note that dispersal occurs in pulses. The dispersal of An, 
arabiensis is slower than that of An. gambiae s.s., prob- 
ably because of the drier conditions in the Sahel and 
An. gambiae s.s. reached the area before An. arabiensis 



(Figure 15). The simulations show that establishment in 
an already occupied area is a much slower process com- 
pared to the case of no competition. From the simulations 
we can also speculate on whether the dominance of one 
species can act as a barrier to genetic flow, like a mountain 
range or dessert. This also raises some questions regarding 
whether hibernation or dispersal is the mechanism behind 
the dominance of the An. gambiae s.s. M form in parts of 
Mali. Although there are strong indications that the An. 
gambiae s.s. M undergoes aestivation during the dry sea- 
son [57,58,71], it is also possible that the persistence of 
the An. gambiae s.s. M form in the Niono district in Mali 
can serve as a refuge during the dry season [153]. In both 
cases the M form receives a kick-start at the beginning 
of the rainy season, and might slow down the dispersal 
of An. arabiensis and the S form of An. gambiae s.s. A 
similar mechanism could contribute to the dominance of 
An. arabiensis in Ethiopia in the Turkana district, where 
the presence of An. arabiensis prevents rapid invasion by 
An. gambiae s.s. 
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Figure 15 Number of months required to reachi a density of 20mosquitoes/km^. Panel 1 (left to right) represents /\n. arabiensis in experiment 
mosqTran(a), panel 2, An. gambiae s.s. with the presence of An. arabiensis (mosqTran(a)), and panel 3, An. gambiae s.s. with no competition 
(mosqTran(b)). The red solid circle and triangle indicate the initial position of An. arabiensis and An. gambiae s.s., respectively. 
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Conclusions 

We developed a model to predict the presence and abun- 
dance of An. arabiensis and An, gambiae s.s. The model is 
age-structured and includes mosquito dispersal 

Sensitivity tests showed that as well as temperature, 
relative humidity and mosquito size are important fac- 
tors in malaria transmission. The result for body size is 
in line with several studies [7,51,54,55,88,154] and thus 
the model captures some of the aspects related to higher 
survival among larger individuals. Note that we have not 
accounted for the higher metabolism in large mosquitoes 
[71], which might reduce survival under warm and dry 
conditions. There are also contrasting results with respect 
to body size and egg production [155]. It is likely that 
there is an optimum size that depends on the environment 
and is a function of temperature and humidity. Currently 
there are few results to back up this statement. However, 
Sanford et al. found significant differences in Anopheles 
gambiae s.s. wing length between Mali and Guinea-Bissau 
[156]. 

We show that relative humidity can be important for 
malaria transmission. Several models have neglected the 
role of (relative) humidity [29,157] and it is true that des- 
iccation might not be a driver of mortality at moderate 
humidity (>70%?). The main argument for leaving out 
this parameter is the corresponding reduction in model 
complexity. As long as rainfall drive the carrying capac- 
ity, mosquito numbers will be restricted at lower humidity 
(no rain), and as a consequence the resulting number of 
mosquitoes can be limited for the wrong reasons, but with 
the correct result. For example, Ermert et al. [28] han- 
dle this deficiency by reducing vector survival during dry 
atmospheric conditions, defined as a function of 10-day 
accumulated rainfall. More studies on the survival of An. 
gambiae s.l. in relation to size and relative humidity in the 
range 5-40% are needed for more confidence in the role of 
humidity in the survival oiAn. gambiae s.l. 

Assumption of exponential mortality has several advan- 
tages (see Figure 5 for examples of models in which 
exponential mortality is used). The model becomes fast to 
solve and it is easier to analyse the equations analytically. 
However, several studies have shown that mortality oiAn. 
gambiae s.l. is not exponential, and that inclusion of an age 
dimension alters the expected outcome of interventions 
targeted to reduce the vector population [50]. Therefore, 
we believe that models in which age-dependent mortal- 
ity is assumed should be further explored. The sensitivity 
tests also suggest that carrying capacity within a restricted 
area plays a role in the distribution of An. arabiensis and 
An. gambiae s.s. The true carrying capacity is hard to esti- 
mate on a continental scale and thus relies on qualified 
guesswork taking into account rainfall, groundwater and 
soil saturation, for example. Carrying capacity influences 
not only the relative distribution of the two species but 



also the total number of mosquitoes. To correctly estimate 
the biting rate, a correct estimate of carrying capacity is 
required, and thus more work is needed to parametrize 
puddle formation. It should also be noted that no current 
large-scale models can describe the formation of puddles 
as rivers retreat, as described by Animut et al. [158]. 

Experiment pBlood2D showed how the model responds 
to the parameter P{B), the probability of finding a blood 
meal. P{B) is important in describing a realistic distri- 
bution of An. arabiensis and An. gambiae s.s. Thus, we 
hypothesize that the large-scale distribution of bovines is 
key to the success of An. arabiensis. Likewise, large-scale 
human density favours the presence of An. gambiae s.s. 

Finally, experiment mosqTran showed how the initial 
conditions influence the dispersal of An. gambiae s.s. 
(and An. arabiensis). The distribution oiAn. gambiae s.s. 
changes dramatically with the presence of An. arabiensis, 
and thus the initial model conditions are highly rele- 
vant for correct description of the distribution of the two 
species. When rainfall is highly seasonal, the first come, 
first served principle seems to be important for the suc- 
cess of a species in drier conditions. Whether or not this 
plays a role in the evolution of aestivation in An. gam- 
biae s.s. M form [57] is a question that should be further 
investigated. 

The strong influence of initial conditions on dispersal of 
the An. gambiae complex is not irrelevant when assess- 
ing the impact of climate change, since vectorial capacity 
varies between species. 

The availability of mosquito models allows researchers 
to build on and improve our understanding of the role 
of the An. gambiae complex in malaria transmission. We 
hope to refine the model as new data on mosquito biol- 
ogy become available, and to incorporate the effects of 
interventions. 
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